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Abstract 

Meta-regression models are commonly used to synthesize and compare effect sizes. 
Unfortunately, traditional meta-regression methods are ill-equipped to handle the com¬ 
plex and often unknown correlations among non-independent effect sizes. Robust variance 
estimation (RVE) is a recently proposed meta-analytic method for dealing with depen¬ 
dent effect sizes. The robumeta package provides functions for performing robust variance 
meta-regression using both large and small sample RVE estimators under various weight¬ 
ing schemes. These methods are distribution free and provide valid point estimates, 
standard errors and hypothesis tests even when the degree and structure of dependence 
between effect sizes is unknown. 


1. Introduction 

Meta-analysis refers to a collection of quantitative methods for combining evidence across 
studies (Hedges and Olkin 1985). Often this evidence is defined in terms of effect sizes, 
which are standardized so as to be comparable across studies. It is reasonable to treat these 
effect sizes as independent when estimators are derived from independent experiments or non¬ 
overlapping subject pools. However, non-independent effect sizes and effect size estimates are 
ubiquitous in published research (Van den Noortgate, Lopez-Lopez, Marin-Martinez, and 
Sanchez-Meca 2013) and ignoring or misspecifying their covariance structure can result in 
estimates that are not valid (Matt and Cook 2009). 

Two common forms of dependence that arise in meta-analysis are correlated estimation errors 
and correlated effect size parameters. Consider the simple model 

Ti = 0i + £j ( 1 ) 

where Ti is an estimate of the i th effect size, 9j is the true effect size, and S{ is the error 
resulting from our estimate of 0* with T). When 6i is fixed (or varies only as a function of 
known characteristics), dependence can occur through the estimation errors (ejs). In practice, 
this form of dependence is common to studies where multiple estimates are obtained from the 
same units (e.g., multiple outcomes on the same individuals, multiple time points), or where 
a common control group is used for multiple experimental contrasts. 

When the 6iS are considered random, the dependence can occur through the 0iS. Dependence 
among effect size parameters is commonly encountered when a hierarchical relationship exists 
among studies, such as reports published by the same laboratory, different groups using the 
same equipment, or multiple experiments reported in a single study. In this case, it is the 
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studies which are nested within a larger category. For a more detailed treatment of effect size 
dependency in meta-analysis see Gleser and Olkin (2009) or Becker (2000). 

The problem with dependence is that traditional meta-analysis models assume independent 
effect sizes. When studies report multiple effect sizes, two approaches are common. One is to 
the use study average effect sizes, which result in a loss of information. A second approach 
is full multivariate meta-analysis (Hedges and Olkin 1985; Raudenbush, Becker, and Kalaian 
1988; Kalaian and Raudenbush 1996). This is the optimal method for handling correlated 
data when the within-study covariance structure is known. Unfortunately, this information 
is rarely reported in published research. Without knowledge of the within-study correlations 
a number of factors should be considered before multivariate methods are used in practice. 
Those interested in a more detailed discussion of this topic should see Jackson, Riley, and 
White (2011). 

Recently, Hedges, Tipton, and Johnson (2010) provided a new approach for handling non- 
independent effect sizes that does not require knowledge of the within-study covariance. The 
RVE method extends work on heteroskedasticity-robust (Eicker 1967; Huber 1967; White 
1980) and clustered (Liang and Zeger 1986) standard errors in the general linear model, to 
the forms of dependence and weighting most often encountered in meta-analysis. The RVE 
approach has a number of desirable qualities. In meta-regression, the coefficients obtained 
from RVE are consistent estimates of the underlying population parameters under a broad set 
of conditions including non-normality. The results do not require the predictor variables to 
be fixed, as is the case in traditional regression methods used in meta-analysis (Hedges 1982). 
Additionally, RVE produces valid standard errors, point estimates, confidence intervals, and 
significance tests when effect sizes are non-independent without needing to model the exact 
nature of this dependence. The RVE approach is increasingly being used in psychology, 
intervention research, medicine and ecology (see Tipton in press). An overview of RVE 
methods geared towards practitioners, including a software tutorial for the Stata program 
(Hedberg 2011) and SPSS, can be found in Tanner-Smith and Tipton (2013). 

A previous limitation of RVE is that its application was limited to situations where the 
number of studies is large; this is because the results described in Hedges et al. (2010) depend 
on asymptotic convergence. Simulation studies on the rates of convergence for regression 
coefficients (Williams 2012), the standardized mean difference (Hedges et al. 2010), as well as 
the log odds ratio, log risk ratio, and risk difference effect sizes (Tipton 2013), suggested that 
while confidence interval coverage for intercepts was adequate with as few as 10 studies, the 
coverage rates for slope intervals were only nominal when the number of studies was greater 
than 40. This is problematic since work by Ahn, Ames, and Myers (2012) and Polanin (2013) 
suggests that more than 50% of meta-analyses in psychology and education have fewer than 
40 studies. 

More recently, Tipton (in press) provided adjustments to the RVE estimators and degrees 
of freedom that greatly improve its performance in small samples. These are particularly 
important when the number of studies is less than 40, but also when the covariates are unbal¬ 
anced or highly skewed. Tipton (in press) shows that it is hard to know a priori when these 
corrections are needed and suggest implementing them in all RVE analyses. The robumeta 
package provides users with the ability to perform meta-regression with robust standard errors 
using both the original large-sample (Hedges et al. 2010) and small-sample adjusted (Tipton 
in press) procedures. We briefly review these developments in RVE before proceeding to an 
example of RVE meta-regression using robumeta. 
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2. RYE 


Suppose we have j = 1,..., m studies, each with with kj effect size estimates, T r Let study 
j also have a vector of kj residuals ej, a design matrix Xj, and weight matrix W j. The linear 
model relating the components can be written as 

T = X/3 + e (2) 

where the vectors T = (T , 1 ,...,T' m )' and e = {e\, ... £' m )' each contain kjxl stacked vectors 
corresponding to our m studies, design matrix X = (X^,..., X^/ contains m stacked kj x p 
matrices, and (3 = {J3\,...,(3 P )' is a p x 1 vector of unknown regression coefficients. For a 
given block diagonal matrix of weights W = diag(W / 1 ,..., W' n )', the weighted least squares 
estimates of (3 are 

( m \ / m \ 

X V " X j { » " 1 j (3) 

When the covariance of the £j in study j is Xj, the exact pxp covariance matrix for b is 


( m 

£ X'w-x. 

The problem that arises in estimating V(b) stems from difficulties in estimating the within- 
study covariance matrix Xj. Although the individual variances are known, calculating the 
covariances of £ :i in Xj depend on correlations that are often unreported. In RVE, the cross 
products of residuals for study j, e^e'-, are used as a rough estimate of the unknown Xj. In 
doing so the estimate of V(b) becomes 


£ 


X'WjXjWjXj 



X'WjXj 


(4) 


yi?. 


\ 
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|vx’w x j 


\ 

E x i w i x iJ 


(5) 


Although the e 3 ef are rather poor estimates of the individual Xj’s, the V R estimator converges 
to the correct value, V(b), as the number of studies in the meta-analysis m —> oo (see Hedges 
et al. 2010, Appendix). Inferences made on the regression coefficients are based on these 
robust standard errors. 


2.1. Weights 

In traditional meta-analysis, inverse variance weights play two roles, First, they produce the 
most efficient estimate of (3. Second, the statistical estimator of the variance requires the 
weights to be inverse; when they are not, the estimate of the variance can be severely biased. 
In RVE, inverse variance weights are also used, but here they are only important for efficiency. 

Weighting improves statistical efficiency by allocating more weight to studies that are more 
precise (i.e., have a smaller variance). Importantly, RVE provides asymptotically accurate 
estimates of the standard errors and valid inferential methods for any set of weights. This 
means the question of weights in RVE is only about efficiency. The most efficient weights would 
be Wj = XJ 1 . Clearly this is problematic since RVE is used when XJ 1 is unknown. For this 
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reason, in RVE weights are developed based on a working model of the unknown covariance 
structure Y, a j so that W j ~ X” 1 . Two working models (hierarchical and correlated effects) 
have been proposed to model the types of dependency most commonly encountered in rneta- 
analytic research. Both of these weighting methods are only approximately efficient and in 
practice, researchers are urged to choose the weighting method based on the most common 
source of dependence in their data. 

Correlated effects 

For study j, a working model based on the correlated effects dependence structure can be 
written as 


Fk/ji — w — *Tjr' T puj (J j Ij) T Uj Ij (6) 

where Jj is a kj xkj all-ones matrix, I j is a k j xkj identity matrix, vj is the estimation error 
variances for the kj effect sizes, and p is a common correlation between the kj effect sizes. Since 
the correlated effects model assumes dependence arises from measurements made on the same 
number of subjects (eQ, it is assumed the error variances Vij would be roughly constant within 
studies = Vj. Therefore, the simplified working model assumes a common correlation 
between within-study effect sizes. Approximate inverse variance weights are calculated by 
first partitioning the total study effect size variance equally among the kj effect sizes, and 
then taking the inverse. More generally, we do not assume the error variances are constant 
within studies i Zy f Vj, and approximate inverse variance weights are calculated using 


Wj = Wjlj = {l/[kj(v.j + r 2 )]}Ij (7) 

where all the effect sizes in study j are assigned the same weight Wij = Wj , v.j is the average 
study variance, and r 2 is an estimate of the between-study variance in study-average effect 
sizes. 

To estimate t 2 , we compute an initial meta-regression where all the effect sizes within study 
j are given equal weights (wj = 1 /(kjVj)). Thus, the weighted residual sum of squares Qe is 


Qe = Et'w/t,- ' ^x'w/x, 


3 =1 


C =1 


C =1 


E X ^ T 

J =1 


3 > 


( 8 ) 


and the method of moments estimate of between studies variance f 2 is given by 


f 2 = 


Q E - m + tr [v (E^r IfX'X, 


J2T=i k 3 w i ~ tr 


+ 


Ptr [V (E” : 


= 1 kn 


X'JjXj - X'X. ; 




(9) 


Note that this depends on an unknown value for the common correlation p. In practice, these 
results do not vary much for p £ (0,1); later we will introduce a sensitivity approach that can 
be implemented in practice. 
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Hierarchical effects 

For study j, a working model based on the hierarchical dependence structure can be written 
as 


51 a j — T T Vj (10) 

where Jj is a k j x k j all-ones matrix, I j is a k j x k j identity matrix, Vj is a kj x kj is a diagonal 
matrix of error variances for study j, co 2 is a measure of effect size variation within study j, 
and t 2 is a measure of effect size variation between studies. Approximate inverse variance 
weights for study j are retrieved from 

Wj = diag(u>ij,..., w k j) (11) 

where Wij = +t 2 +<u 2 ). The weights are only constant in study j if the error variances 

for each effect size are equal. 

To estimate the parameters t 2 and w 2 , we first obtain the weighted residual sum of squares us¬ 
ing (8). We also calculate an additional sum of squares using the results from our preliminary 
meta-regression 


Qi = E ( T i - x ; - pj' J i (t? - x . - p) • ( 12 ) 

3 = 1 

Hedges et al. (2010) provide the following method of moments estimators for r 2 and cj 2 


A 2 (Qi-C 1 )-A 1 (Q e -C 2 ) 

B1A2 — B2A1 


(13) 


_ Q E — C2 ^ B2 

T2 U 2 T 2 

Details on these equations can be found in Hedges et al. (2010). 


(14) 


3. Small-sample adjustments 

As previously mentioned, while e 3 e( is a crude estimator of Hj for individual studies, 
converges in probability to V(b) asymptotically (Hedges et al. 2010). The properties, and 
limitations, of linearization estimators like V 2? in small samples are well documented, and a 
number of solutions have been proposed to address these limitations. Interested readers can 
see Irnbens and Kolesar (2012) for a discussion of these solutions. 

In this section, we review two types of small sample adjustments proposed by Hedges et al. 
(2010) and extended by Tipton (in press). The first of these corrections is to the RVE 
estimator itself, while the second is to the degrees of freedom used for making inferences 
regarding the meta-regression coefficients (3. 
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3.1. Covariance structure adjustments 

The first adjustment is to the estimator which can be written more generally as 

X'W ; X ; j ^ Xj- Wj AjOjc'j Aj Wj X j 

where in general A jej = Aj(I — H )je, H = XQX'W, Q = (X'WX) -1 , and e is a £k j x 1 
vector of true residuals. Hedges et al. (2010) proposed that this estimator could be improved 
by using 



— 1 


x'w ; x, 


( 15 ) 



A f TJ = [ m / ( m — p)] 1 / 2 ^- (16) 

where I j is a kj by kj identity matrix. Simulations indicate that this adjustment is inadequate 
except in very specific cases and Tipton (in press) develop a new solution based on extensions 
to the Bias Reduced Linearization (BRL) method proposed by McCaffrey, Robert M. Bell, 
and Botts (2001). Here since V* also depends on W and working covariance matrices S aj -, 
the Skj by Sk ? - adjustment matrix A j is specified separately for the hierarchical 


A f = w: 1 / 2 [w: 1/2 (i - H i 3 )w : 3/2 ]- 1/2 W : 1/2 

and correlated effects weighting model 


(17) 


Af = (I-H jj)- 1 /*, (18) 

where Hjj = XjQX'-Wj. Additionally, if non-efficient weights are required (in robumeta 
non-efficient weights can be specified using the userweights option) the adjustment matrix 
is specified as follows 


Af = ,.f 


(I-H).V(I-H) 


,1-V 2 


(19) 


where (I — H)j includes the kj rows of (I^H) associated with study j and the average variance 
of the study j’s effect sizes, u.j = (1 /kj) Yivij. 

In the above, for a general matrix B note that I = B^BB' 1 / 2 , and B -1 / 2 = PA" 1 / 2 ?' 
where PAP' is the eigen decomposition. See the Appendix of Tipton (in press) for a detailed 
account of the development of these adjustments and their use with non-inverse variance 
weighting schemes. 


3.2. Degrees of freedom adjustment 

The second adjustment proposed by Hedges et al. (2010) was to use the t-distribution with 
df = m — p when making statistical inferences using . Building off of Bell and McCaffrey 
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(2002), however, Tipton (in press) suggests using the Satterthwaite approximation (Satterth- 
waite 1946) to estimate the distribution of V* with Xdf Sk i where the degrees of freedom dfsk 
is equal 2/c^ 2 , and cv is the coefficient of variation. 

These degrees of freedom are calculated separately for each regression coefficient under study¬ 
ing using 


fe = £^) 2 /E 4 (20) 

where A jk are the eigenvalues of W~ 1//2 (^() gj/tgj/jW -1 / 2 where g jk = (I — H)' AjWjXjQ4- 
Here Ik is a vector of 0’s and l’s that indicates the regression coefficient, A j and W j are 
defined as before. Further details can be found in Tipton (in press). Simulations by Tipton 
(in press) show that the Satterthwaite approximation is valid so long as the df > 4. This 
means that RVE can be effectively used to make inferences even with a very small number of 
studies. 


4. Testing and confidence intervals 

Robust tests of the null hypothesis Ho : /3& = 0 are conducted using the test statistic 

ti = h/VWk (2i) 

where bk is the estimate of /?&, is the k th diagonal of the V* matrix, and where we reject 
Ho : /3fc = 0 when |t^| > tdf, Q - Additionally, a 95% confidence interval for the kth regression 
coefficients /5k is given by 


bk — tdf,ct^Vlk — Pk < bk + tdf,a\/Vkk- ( 22 ) 

Under the original Hedges et al. (2010) approach, V£ fc is calculated using equation 5 and 
df = m — p. When the small sample corrections from Tipton (in press) are used, is 
calculated using A^, A*p, A) v for the correlated effects, hierarchical effects and non-efficient 
weighting models, respectively. Here the degrees of freedom dfsk , as specified in equation 20, 
vary from covariate to covariate. 


5. The robumeta package 

The robu () function is the main fitting function for robumeta and implements rneta-regression 
using the large and small-sample RVE estimators described in Hedges et al. (2010) and Tipton 
(in press), respectively. The arguments available in robu() are 

robu(formula, data, studynum, var.eff.size, userweights, modelweights = 
c("C0RR", "HIER"), rho = 0.8, small = TRUE,...) 

The formula object is similar to that found in other linear model specification routines (e.g., 
lm(), glm()). A typical formula will look similar to y ~ xl + x2. . ., where y is a vector 
of effect sizes and xl + x2. . . are the user-specified covariates. An intercept only model 
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can be specified with y ~ 1 and alternatively, the intercept can be omitted by specifying y 
-1 and multiplicative interaction terms can be specified with xl*x2. . .. The data argument 
indicates the user-specified data frame which contains the following 3 columns 

1. Effect sizes (effect. size): User-specified column of effect sizes calculated from original 
research reports. 

2. Effect size variances (var. ef f . size): User-specified column of within-study variances 
for each effect size. 

3. Study identification (studynum): User-specified column of study unique values used to 
identify studies. 

4. User-specified weights (userweights): User-specified weights if non-efficient weights are 
of interest (optional). 

modelweights is the user-specified option for selecting either the hierarchical (HIER) or cor¬ 
related (CORR) effects model. By default the correlated effects model is used, rho is the 
user-specified within-study effect size correlation, which must be specified when using the 
correlated effects model. The default value for rho is 0.8 and user-specified values must be 
between 0 and 1. small = TRUE or small = FALSE allow the user to specify whether or not 
the Tipton (in press) adjusted estimation procedures will be used to fit the model, small = 
TRUE is the default. The robu() function returns a fitted object of class robu which contains 
a list of summary information calculated during the model fitting procedure. The print () 
method provides a formatted output containing this summary information. 


6. Examples 


6.1. Example 1 

In the first example we use data from a meta-analysis conducted by Oswald, Mitchell, Blanton, 
Jaccard, and Tetlock (2013) examining the predictive validity of the Implicit Association Test 
(IAT) and various explicit measures of bias for a variety of criterion measures of discrimination. 
The oswald2013 dataset provided in robumeta contains 308 effect sizes collected from m = 46 
individual studies. 

For this example, we will look at a subset of the oswald2013. exl dataset which contains only 
the studies where the criterion measure of discrimination used was either neurological activity 
or response latency. The (oswald2013. exl) dataset contains 32 effect sizes collected from 
m = 9 individual studies. Additionally, the reported correlations have been transformed to 
Fisher’s Z r (Gleser and Olkin 2009). As the majority of effect-sizes in our dataset result from 
each study including multiple measures taken on the same individuals, we will use correlated 
effects weights to fit our model. 

6.2. Intercept-only model 

When performing a meta-regression it is standard practice to first calculate an average effect 
size without conditioning on the study covariates. In RYE this can be followed with a sen- 
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sitivity analysis which will be discussed in more detail in the following section. The fitting 
function for the intercept-only model is specified with the robuQ function 

oswald_intercept <- robu(formula = effeet.size ~ 1, data = oswald2013.exl, 

studynum = Study, var.eff.size = var.eff.size, 
rho = .8, small = TRUE) 

and the results displayed with print () as follows 

print(oswald_intercept) 

RVE: Correlated Effects Model with Small-Sample Corrections 
Model: effect.size ~ 1 

Number of studies = 9 

Number of outcomes = 32 (min = 1 , mean = 3.56 , median = 2 , max = 11 ) 

Rho = 0.8 

I.sq = 82.05272 

Tau.Sq = 0.1849812 

Estimate StdErr t-value df P(|t|>) 95% CI.L 95% CI.U Sig 
1 intercept 0.277 0.181 1.53 7.84 0.164 -0.141 0.695 

Signif. codes: < .01 *** < .05 ** < .10 * 

Note: If df <4, do not trust the results 

In addition to the coefficient table, the output from print () contains the formula call used 
to fit the model, summary information for the total number of studies and outcomes included 
in the analysis, and a coefficient table. In addition, output for the correlated effects model 
contains the user-specified p value, the estimate of between-study variance in study-average 
effect sizes r 2 , and a descriptive statistic for the ratio of true heterogeneity to total variance 
across the observed effect sizes, / 2 (Higgins, Thompson, Deeks, and Altman 2003). 1 2 is 
calculated using 


I 2 = [(Qe - rf/)/Q E ]100% (23) 

where Qe is defined as the weighted residual sum of squares as in equation 8. Unlike r 2 , 1 2 
is not on the same scale as the effect size index used in the analysis and is better suited to 
anwsering questions regarding the proportion of observed variance resulting from real effect 
size differences. 

6.3. Sensitivity analysis 

As mentioned previously, an estimate of r 2 is required to develop approximately efficient 
weights for the meta-regression model. The method of moments estimator for r 2 requires a 
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value is specified for p. For this reason Hedges et al. (2010) suggests conducting a sensitivity 
analysis to determine the effect of p on r 2 . Generally, results from simulation studies suggest 
that t 2 and the meta-regression coefficients are relatively insensitive to changes in p (Ishak, 
Platt, Joseph, and Hanley 2008; Williams 2012; Tipton 2013); however, if between-study co- 
variances were to be considerably smaller than within-study covariances, these results might 
not hold. It should also be noted that a conservative approach and an external information 
approach have also been advocated for estimating p in terms of developing efficient weights. 
The conservative approach involves setting p to 1, which ensures studies don’t receive addi¬ 
tional weight due to having more effect sizes. Additional information on these strategies can 
be found in Hedges et al. (2010). 

The sensitivtyO function in robumeta computes r 2 , the average effect size and associated 
standard error for different values of p in the interval (0,1). 

R > sensitivity(oswald_intercept) 



Type 

Variable 

& 

o 

II 

o 

rho=0.2 

o 

II 

o 

5h 

rho=0.6 

& 

o 

II 

o 

00 

rho=l 

1 

Estimate 

intercept 

0.277 

0.277 

0.277 

0.277 

0.277 

0.277 

2 

Std. Err. 

intercept 

0.181 

0.181 

0.181 

0.181 

0.181 

0.181 

3 

Tau.Sq 


0.184 

0.184 

0.184 

0.185 

0.185 

0.185 


From the sensitivity analysis it seems t 2 and subsequently the average effect size are relatively 
robust to different values of p. 

6.4. Forest plot 

In meta-analysis, forest plots provide a graphical depiction of effect size estimates and their 
corresponding confidence intervals. The f orest. robu() function in robumeta can be used 
to produce forest plots for RVE meta-analyses. The function requires the grid (R Core 
Team 2013) package and is based on examples provided in Murrell (2011). As is the case in 
traditional forest plots, point estimates of individual effect sizes are plotted as boxes with areas 
proportional to the weight assigned to that effect size. Importantly, here the weight is not 
necessarily proportional to the effect size variance or confidence intervals, since the combined 
study weight is divided evenly across the study effect sizes. Two-sided 95% confidence intervals 
are calculated for each effect size using a standard normal distribution and plotted along with 
each block. The overall effect is included at the bottom of the plot as a diamond with width 
equivalent to the confidence interval for the estimated effect. 

The RVE forest function is designed to provide users with forest plots which display each 
individual effect size used in the meta-analysis, while taking into account the study- or cluster- 
level properties inherent to the RVE analysis. As such, the user must specify columns from 
their original dataset that contain labels for the study or cluster and for the individual effect 
sizes. For instance, labels for the study column might be author names with corresponding 
publication years, and would be specified in the forest.robu function using the study.lab 
argument. Labels for individual effect sizes might be “Math Score” or “Reading Score” for a 
meta-analysis that included such measures or as simple as “Effect Size 1” and “Effect Size 2,” 
and can be specified with the es. lab argument. Any number of additional columns can be 
specified to be plotted along side the confidence interval column and can be specified with 
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the following syntax "argl" = "arg2" where "argl" is the title of the column on the forest 
plot, and "arg2" is the name of the column from the original data.frame that contains the 
information to be displayed. An RVE forest plot with study labels corresponding to the 
author name and publication year, and effect size labels corresponding to the criterion for 
the Oswald data (see Figure 1 on page 12) is specified using the syntax below. In addition, a 
column of the RVE weights titled “Weight” and a column of effect sizes titled “Effect Sizes” is 
also included to illustrate the syntax for adding additional columns. 

forest.robu(oswald_intercept, es.lab = "Crit.Cat", study.lab = "Study", 

"Effect Size" = effeet.size, # optional column 
"Weight" = r.weights) # optional column 


6.5. Example 2 

In the following section we walk through an example where robumeta is used to fit simulated 
data described in Tanner-Smith and Tipton (2013). The dataset in question resulted from 
a fictional meta-analysis that examined the effectiveness of a brief alcohol intervention on 
alcohol intake. Intervention effectiveness was determined by differences in alcohol consump¬ 
tion between treatment and controls following the intervention. In this example standardized 
mean difference effect sizes are provided using Hedges’ g, however, any effect size can be used 
with the robu() function. The hierdat dataset contains m = 15 unique treatment centers 
with 68 effect sizes. The covariates of interest in this example are followup, which indicates 
the length of the follow-up period before alcohol consumption was measured, in months and 
binge which indicates whether or not the alcohol consumption measure targeted binge drink¬ 
ing specifically. As the majority of dependency in this meta-analysis resulted from studies 
nested within treatment centers, we use the hierarchical effects weights to fit our model. 

6.6. Covariate levels 

In RVE dependencies result from hierarchically structured data. Borrowing terminology from 
the multilevel and hierarchical linear modeling literature, this hierarchical structure takes 
the form of lower-level observations nested within higher-level groupings. Within the RVE 
framework we are generally concerned with two-level structures where level-1 (the lowest 
level) represents effect sizes in the correlated effects case and studies in the hierarchical effects 
case, while level-2 represents the study-level in the correlated effects case and a higher-order 
grouping level such as research lab in the hierarchical effects case. 

Since covariates can vary both between-level-2 and within-level-2 it is an important pre¬ 
processing step in any RVE analysis to identify the nature of this variation prior to model 
fitting. For example, at first glance the followup variable appears to vary both within- and 
between-level-2. This indicates the length of follow-up periods varies both across studies 
conducted by the same treatment center, and then between the treatment centers themselves. 
Failure to distinguish between these effects can complicate interpretation later on as the 
regression coefficient actually represents a weighted combination of both between- and within- 
level-2 effects, which can be problematic when the effects differ in direction or magnitude. 
For instance, if treatment centers which produce higher quality studies also tend to employ 
longer follow-up periods the effect of follow-up length on treatment effect sizes might differ 
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Forest Plot 

RVE: Correlated Effects Model 
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Figure 1: This figure shows a forest plot for the Correlated Effects model. 
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when examined within a given treatment center’s studies, compared to the effect of follow-up 
length across treatment centers. 

To partition these two types of effect we can create group mean X,j and group centered Xfj 
versions of followup using the convenience functions group.mean() and group, center () 
provided in robumeta. 

R > hierdat$followup_m <- group.mean(hierdat$followup, hierdat$studyid) 

R > hierdat$followup_c <- group.center(hierdat$followup, hierdat$studyid) 

Subsequently, including both followup_m and followup_c in our meta-regression model as 
follows 


Tij =/3 0 + X mj P! + Xgfo + ■■■ (24) 

allows us to decompose the effects of follow-up length into two independent estimates where 
B\ and B 2 represent the effect of a 1-month increase in A,,- and a 1-month increase in Xfj 
on Tij , respectively. Another important feature of centering is that it often indicates that a 
variable only varies at the study level, not effect size level. We therefore encourage researchers 
to examine these variables before proceeding with their analysis. 

6.7. Meta-regression model 

Now that we have completed the necessary data pre-processing we can fit our rneta-regression 
model using hierarchical effects weighting with the following command 

R > model_hier <- robu(effectsize ~ followup_c + followup_m + binge, 

data = hierdat, modelweights = "HIER", studynum = 
studyid, var.eff.size = var, small = TRUE) 

The print () method can be used to present the results 

R > print(model_hier) 


RVE: Hierarchical Effects Model with Small-Sample Corrections 
Model: effectsize ~ followup_c + followup_m + binge 


Number of clusters = 15 

Number of outcomes = 68 (min = 1 , mean =4.53 , median = 2 , max = 29 ) 
□mega.sq = 0.1650524 
Tau.Sq = 0.02479249 


Estimate StdErr t-value df P(|t|>) 95% CI.L 

1 intercept -0.154226 0.147671 -1.044 6.07 0.3361 -0.51461 

2 followup_c -0.000162 0.000695 -0.233 1.30 0.8470 -0.00534 

3 followup_m 0.003467 0.002320 1.495 3.28 0.2243 -0.00357 


95% CI.U Sig 
0.20615 
0.00502 
0.01050 
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4 binge 0.666645 0.115623 5.766 4.33 0.0035 0.35514 0.97815 *** 

Signif. codes: < .01 *** < .05 ** < .10 * 

Note: If df <4, do not trust the results 


The results of the meta-regression suggest the p-value for followup_c and followup_m are 
not reliable, as the Satterthwaite degrees of freedom are less than 4. This may result from the 
fact that center is contributing more than half of the effect sizes included in our meta-analysis 
did not actually vary the followup interval across studies, resulting in an unbalanced covariate. 
This is in line with results from simulation studies which suggest the Satterthwaite degrees 
of freedom are typically smaller for unbalanced covariates (Tipton in press). The effect of 
whether or not the study used a binge drinking measure (binge), however, was found to be 
positive and significant at a = 0.05. This suggests that studies which used binge drinking 
as their criterion measure were associated with increased treatment effect sizes compared to 
those studies that used some other measure of alcohol consumption. 


7 . Conclusions 

Robust variance estimation (RVE) is a user friendly procedure grounded in work on 
heteroskedasticity-robust (Eicker 1967; Huber 1967; White 1980) and clustered standard er¬ 
rors (Liang and Zeger 1986). RVE holds promise for alleviating the problems introduced 
by unknown within-study correlations in met a-analysis (Jackson et al. 2011). This potential 
utility of RVE methods has been further expanded recently through the small- sample ad¬ 
justments developed by (Tipton in press) which greatly improve the finite-sample properties 
of RVE. The robumeta package provides functions for conducting meta-regression with both 
the large-sample (Hedges et al. 2010) and adjusted (Tipton in press) RVE estimators and 
provides functions for producing forest plots consistent with the RVE methodology. 
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